# Task to perform the calculation to get the signal-to-noise spectrum per pixel
# and other complementary informations to help the user set up the observation.
# Written by Sergio Scarano Jr (scarano@astro.iag.usp.br) - jul 2009.

import numpy
import scipy

#execfile(os.path.join(_path,"config.py"))
execfile(os.path.join(_path,"shiftspec.py"))
execfile(os.path.join(_path,"configgrating.py"))
execfile(os.path.join(_path,"configreadout.py"))

obs=open(obsfile)
teo=open(teofile)

waveobs=[]
fluxobs=[]

for ss in obs.readlines():
    waveobs.append(float(ss.split()[0]))
    fluxobs.append(float(ss.split()[1].strip()))

waveobs=scipy.array(waveobs)
fluxobs=scipy.array(fluxobs)

waveteo=[]
fluxteo=[]

for ss in teo.readlines():
    waveteo.append(float(ss.split()[0]))
    fluxteo.append(float(ss.split()[1].strip()))

waveteo=scipy.array(waveteo)
fluxteo=scipy.array(fluxteo)/1E16

# Using a spline to interpolate the observational data obtained at the official SOAR
# Telescope site. Since the data is limited to values lower than 8891 A, the remnant
# data, associated to larger wavelengths, were obtained by extending a linear fitting
# performed over the observational data larger than 7530 A.

from scipy.interpolate import splrep, splev
tck = splrep(waveteo, fluxteo)
fluxint = splev(waveobs, tck)

# Defining the relative gain

gainrel=fluxint/(fluxobs)  # ATENCAO Correcao 1/4 deve ser inserida
                             # ======= por binagem.

#pylab.plot(waveobs,gainrel)
#pylab.show()

execfile(os.path.join(_path,"cutspec.py"))

from scipy.interpolate import splrep, splev
tck = splrep(waveobs, gainrel)
try:
    gainsub = splev(cutwave, tck)
except ValueError:
    print '''<p style="width:500px;"> <span style="font-weight:bold;"><span style="color:red">ERROR:</span></span> Your input configuration has resulted in a non valid input array. 
It is possible that the template spectrum has no flux at the specified wavelenght range
or something alike. For instance "Quasar" has little emission in the red part of the spectrum
which prevent us from providing a S/N estimation. If you are really sure that this should result 
in a valid spectrum, please contact the SOAR staff with information about your program. We will 
do our best to help!
</p>'''
    exit(0)
fluxadu=cutflux/gainsub
waveadu=cutwave

# Creating a full pixel resolution Goodman wavelenght array inside the limits
# of the observed spectrum

wavefulres=scipy.linspace(lambdainf,lambdasup,4142)

from scipy.interpolate import splrep, splev
tck = splrep(waveadu, fluxadu)
fluxfulres=splev(wavefulres,tck)

# Taking into account the spectral dispersions for the chosen grating
# to evaluate the resolution element.

execfile(os.path.join(_path,"configgrating.py"))
execfile(os.path.join(_path,"configslit.py"))

elresini=lambda0-slitpix/2*disp
elresfin=lambda0+slitpix/2*disp

# Counting how many ADUs are inside the resolution element around a
# central wavelength previously indicated

ind=numpy.where((wavefulres < elresfin) & (wavefulres > elresini))

wavesn=wavefulres[ind]
fluxsn=fluxfulres[ind]

razgain=gain/gainobs

fluxsne=gain*razgain*fluxsn
wavesne=wavesn

signal=fluxsne.sum()  # ATENCAO Correcao 1/4 deve ser inserida
                      # ======= por binagem.


#import pylab as py
print '''<h1>Result:</h1>

<table border="1">

'''
if snrintgr == 's':
    exptime=(0.5*sndesire+0.5*numpy.sqrt(sndesire**2+4*rdnsobs**2))*sndesire/(gain*razgain*fluxsn.sum()/(1./binning))
    print '<tr><td> Exptime = ', exptime, '</td></tr>\n'
    print '</table>'
else:
     #binning=1.0
     snr=(gain*razgain*fluxsn.sum()*exptime/binning)/numpy.sqrt(gain*razgain*fluxsn.sum()*exptime/binning+rdnsobs**2)
     print '<tr><td> S/R = %.2f</td></tr>\n'%(snr)
	 
     snrvec = fluxfulres# / fluxsn.sum()
     snrvec = (gain*razgain*snrvec*exptime/binning)/numpy.sqrt(gain*razgain*snrvec*exptime/binning+rdnsobs**2)
     mask = snrvec[ind] > snr*0.1
     #print snrvec[ind][mask].mean()
     mask2 = snrvec > 0.0 # snr*0.1
     py.subplot(121)
     py.plot(wavefulres[mask2],snrvec[mask2]/snrvec[ind][mask].mean()*snr)
     py.plot(wavesn.mean(),snr,'ro') 
     py.ylabel('Signal-to-Noise per resolution element')
     py.xlabel('Wavelength [$\\AA$]')

     py.subplot(122)
     py.plot(wavefulres[mask2],snrvec[mask2])
     py.plot(wavesn.mean(),snr,'ro') 
     py.ylabel('Signal-to-Noise per pixel')
     py.xlabel('Wavelength [$\\AA$]')

     py.savefig(os.path.join(_tmp,'tmp01.png'))
#     py.savefig('/home/user/soarcalc/www_lna/tmpfiles/tmp01.png')
     
     print '''</table>
<img src="http://www.lna.br/~soarcalc/tmpfiles/tmp01.png">
     
</body>
</html>''' 
